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Three-dimensional integer partitions provide a convenient representation of codimension-one 
three-dimensional random rhombus tilings. Calculating the entropy for such a model is a notoriously 
difficult problem. We apply transition matrix Monte Carlo simulations to evaluate their entropy 
with high precision. We consider both free- and fixed-boundary tilings. Our results suggest that 
the ratio of free- and fixed-boundary entropies is a f ree / 'o fixed = 3/2, and can be interpreted as the 
ratio of the volumes of two simple, nested, polyhedra. This finding supports a conjecture by Linde, 
Moore and Nordahl concerning the "arctic octahedron phenomenon" in three-dimensional random 
tilings. 
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I. INTRODUCTION 

Since the discovery of quasicrystals in 1984 1 , quasiperiodic tilings (such as Penrose tilings) and random rhombus 
tilings 2,3 have been extensively studied as paradigms of quasicrystal structure. Quasicrystals are metallic alloys 
exhibiting exotic symmetries that are forbidden by usual crystallographic rules: octagonal, decagonal, dodecagonal 
or icosahedral symmetries. Nonetheless, the existence of sharp Bragg peaks in their diffraction patterns demonstrates 
a long-range translational and orientational order. As compared to perfect quasiperiodic tilings, specific degrees of 
freedom called "phason flips" are active in random tilings. Despite their random character, the latter still display the 
required long-range quasiperiodic structure. 

When tiles are appropriately decorated with atoms, random tilings become excellent candidates for modeling 
real quasicrystalline materials 4 . Therefore the statistical mechanics of random tilings is of fundamental interest 
for quasicrystal science. In particular, the configurational entropy of their phason fluctuations contributes to the free 
energy of quasicrystal, and might be a key ingredient in order to understand quasicrystal stability 5 . 

In parallel, random tilings have become an active field of research in discrete mathematics and computer science 
(see 6,7 for instance), and many challenging questions remain open for investigation in this field. 

The relation between random tilings and integer partitions provides an important tool for the calculation of random 
tilings entropy 8-12 . Integer partitions are arrays of integers, together with suitable inequalities between these integers. 
One-to-one correspondences can be established between integer partitions and tilings of rhombi filling specified poly- 
hedra. However, it seems that such strictly controlled "fixed" boundary conditions inflict a non-trivial macroscopic 
effect on random tilings 8,13 , even in the thermodynamic limit, lowering the entropy per tile below the entropy with 
free or periodic boundary conditions. This effect a priori makes difficult a calculation of free-boundary entropies via 
the partition method. 

This boundary sensitivity is well described, for the simple case of hexagonal tilings 6,14 , in terms of a spectacular 
effect known as the "arctic circle phenomenon" : the constraint imposed by the boundary effectively freezes macroscopic 
regions near the boundary, where the tiling is periodic and has a vanishing entropy density. Outside these "frozen" 
regions the entropy density is finite and we call the tiling "unfrozen" . The boundary of the unfrozen region appears 
to be a perfect circle inscribed in the hexagonal boundary. The entropy density varies smoothly within the unfrozen 
region, reaching a maximum equal to the free boundary entropy density at the center. 



This quantitative result has never been generalized to higher dimension or codimension tilings, because its gener- 
alization requires the knowledge of the free boundary entropy if one wants to use the same proof as in the hexagonal 
case. However, thanks to numerical simulations, it has recently been conjectured by Linde, Moore, Nordhal 15 that 
in dimensions higher than 2, the corresponding arctic region should be a polytope itself, with flat boundaries. It 
was further conjectured (by Destainville and Mosseri 16 and independently by Propp 17 ) that in this case the entropy 
density should be spatially uniform and maximal in the unfrozen region. These conjectures renew the interest for the 
partition method since the relation between both entropies becomes amazingly simple in this case. 

On the other hand, except an early Ansatz 9 and some exact numerical results for small tilings 11 (up to about 300 
tiles), almost nothing is known about the entropy of codimension-one tilings of dimension larger than 2. Note also 
that some numerical Monte Carlo simulations provided an estimate of the entropy of (codimension 3) 3-dimensional 
random tilings with icosahedral symmetry 18 . 

The present paper is devoted to a numerical investigation of codimension-one three-dimensional tilings. Thanks to 
a powerful transition matrix Monte Carlo algorithm, we achieve precise estimates of both fixed- and free-boundary 
entropies. The latter is calculated via a modified partition method, which produces tilings with fixed boundaries that 
do not impose any strain to the tilings, thus generalizing a former two-dimensional approach 14 . Comparing both 
entropies, we support the above conjecture with good confidence. 

The paper is organized as follows: section II reviews the relation between random tilings and integer partitions, 
introduces the different boundary conditions considered in this paper, and describes the arctic region phenomenon 
conjecture. Section III describes our Monte Carlo method and our numerical results. Discussions and conclusions are 
displayed in the last section. 



In this paper we consider three-dimensional tilings of rhombohedra which tile a region of Euclidean space without 
gaps or overlaps. A standard method 2 ' 19 to generate tilings of rhombohedra (or of rhombi in two dimensions) consists 
of selecting sites and tiles in a D-dimcnsional cubic lattice, then projecting them into a d-dimcnsional subspace with 
D > d. The difference D — d is called the codimension of the tilings. The class of symmetry of a tiling is determined 
by both its dimension and its codimension and we denote tilings with such a symmetry by D — > d tilings. We consider 
in this paper codimension-one, three-dimensional random tilings (i.e. D = 4 and d = 3). These 4^3 tilings are 
composed of four different rhombohedra. All rhombohedra have identical shapes but they occur in four possible 
different orientations. Each rhombohedron is the projections of one of the four different three-dimensional faces of a 
four-dimensional hypercube. The interested reader can refer to [ n ] for a review on codimension-one tilings. 

The correspondence between codimension-one tilings and solid partitions, a reformulation of the cut-and-project 
method, is analyzed in detail in reference [ n ], and generalized to higher codimensions in reference [ ]. 

We now define three-dimensional solid partitions. Consider a three-dimensional array of sides fci x fc 2 x fc 3 . Fix an 
integer p > 0, called the height of the partition problem. Put non-negative integers in the array, no larger than p, 
with the constraint that these integers decrease in each of the three directions of space. More precisely, if i\, %i and 
^3 are indices attached to the boxes of the array (1 < i a < we denote by tii 1 ^ 2 ^ i the integral variables attached 
to these boxes (the parts). Our partition array contains N p — kik 2 k 3 parts. The partition constraint is 



II. PARTITIONS, TILINGS AND BOUNDARY CONDITIONS 
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FIG. 1. One-to-one correspondence between plane partitions (left), hexagonal tilings (center), and de Bruijn lines (right) in 
two dimensions 9,11 ; De Bruijn lines of two families (among three) have been represented. 

Figure 1 displays the tiling/partition correspondence in two dimensions for easy visualization. Plane partitions 
are k x I arrays of non-negative integers smaller than a height p and decreasing in each row and column (left). This 
correspondence is constructed as follows: stack 3-dimcnsional cubes above the boxes of the partition so that the height 
of each stack (center) equals the value of the corresponding partition box (left). Then project this stacking along the 
(1,1,1) direction of the cubic lattice. Faces of the cubes project to rhombi. The so-obtained rhombus tiling fills a 
hexagon of sides k, I and p (center). 

Following the same construction in three dimensions, four-dimensional hypercubes are stacked above the three- 
dimensional partition array, with the heights of the stacks equal to the corresponding parts. Then project into three 
dimensions along the (1,1,1,1) direction of the hypercubic lattice. Like in the hexagonal case, the so-obtained tilings 
fill a polyhedron, a "rhombic dodecahedron" (RD) of integral sides k\, k 2 , k 3 and p (see the outer frame in figure 3). 
The total number of tiles, 

N t = hk 2 k 3 + k x k 2 p + k x k 3 p + k 2 k 3 p. (3) 

We call tilings with rhombic dodecahedron boundaries "i?Z? £?-tilings" and denote their configurational entropy per 
tile 20 by a fixed . 

The source of configurational entropy can be easily understood in either the tiling or the partition representation. 
Figure 2 illustrates the basic "tile flip" in both the 3 — > 2 and the 4 — > 3 tilings. In each case interior tiles are 
rearranged without disturbing the surface of a region. In terms of the equivalent partition, the height of one partition 
element increases or decreases by one unit. This elementary local move is ergodic, and every legal tiling can be reached 
by a succession of such flips. 




(a) (b) 

FIG. 2. Examples of single vertex flips: (a) rotation of 3 rhombi inside hexagon in 3 — > 2 tiling; (b) rotation of 4 rhombohedra 
inside rhombic dodecahedron in 4 — > 3 tiling. Each rotation increases or decreases partition height by 1 unit, as shown. 

An alternative description of the same tilings and partitions is in terms of de Bruijn lines 21 and membranes. Figure 1 
(right) shows the de Bruijn lines of the tiling (center). These are formed by connecting centers of parallel edges on 
every rhombus. There arc three families of de Bruijn lines, each family with an average orientation perpendicular to 
the tile edges. Although de Bruijn lines within a family never cross, de Bruijn line crossings among different families 
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occur at the center of each rhombus. For three-dimensional tilings, we have instead de Bruijn membranes. Elementary 
flips such as those illustrated in Fig. 2 bend the membrane locally, and membrane fluctuations become the source of 
configurational entropy. 

Finally, we note that every tiling can be represented as a d-dimensional directed hypersurface embedded in a D- 
dimensional space. Beginning with d = 2 tilings, this hypersurface consists of the faces of the stacking of cubes visible 
from the (1,1,1) direction (figure 1 (center)). When projected along the same direction on the "real" space E, it 
becomes a plane tiling. It is "directed" because no tile overlaps occur during the projection. The same relationship 
holds for d = 3 tilings: the directed hypersurface consists of the three-dimensional faces of the stacking of hypercubes 
viewed along the (1,1,1,1) direction. Since the hypersurface is directed, it can be seen as a single-valued, continuous, 
piecewise linear function ip from the real space E = R 3 to R. The value of <j> is just the height along the (1,1,1,1) 
axis. In the thermodynamic limit, these piecewise linear functions ip can be coarse-grained to obtain smooth functions 
<f> : R 3 — > R which only contain large-scale fluctuations of the original corrugated membranes 3,14 . 



B. Boundary conditions 

Polyhedral boundary conditions, such as the rhombic dodecahedron bounding RDB tilings, have macroscopic effects 
on random tilings. In the "thermodynamic limit" of large system size, the statistic ensemble is dominated by tilings 
which arc fully random only inside a finite fraction of the full volume and are frozen in macroscopic domains. By 
frozen, we mean they exhibit simple periodic tilings in these domains with a vanishing contribution to the entropy. 
In two dimensions, this is known as the "arctic circle phenomenon", as described in introduction 6,14 . 

Such boundary conditions arc not very physical. Indeed, even if one can imagine situations where the quasicrystal 
is constrained by a flat interface {e.g. growth experiments on a crystalline substrate), the previous considerations 
rely on the assumption that, in the physical quasicrystalline material, the tiles are elementary unbreakable structures. 
However, the system could possibly lower its free energy by breaking some tiles (a line of tiles in 2D or a surface in 
3D), in order to adapt to the constraint, and thereby free the remainder the tiling from the constraint. For sufficiently 
large tilings a net lowering of free energy results. For example, the energy cost of a broken line of tiles grows linearly 
with the length of this line, whereas the free energy difference between a constrained tiling and a free one grows like 
the number of tiles and therefore like length squared (in two dimensions). In the thermodynamic limit, even if it costs 
a great amount of energy to break a tile, the system will eventually prefer to pay this cost. Therefore, the "true" 
thermodynamic entropy density is the free-boundary entropy 22 . 

Consequently, it is desirable to relate fixed boundary condition entropies to the more physical free boundary ones. 
Fortunately, there exists an exact formal relation between these entropies . A quantitative relation has even been 
calculated in the hexagonal case 6,14 , but it has not been possible (so far) to extend quantitative relationships to more 
general classes of tilings. However, a conjecture by Linde, Moore and Nordahl 15 has recently brought some new hope 
(see section II C). 

To exploit the calculational advantages of a partition representation, while achieving the physical free-boundary 
entropy in the thermodynamic limit, we adapt the partition method so that the corresponding tilings exhibit no frozen 
regions. The new boundary, even though fixed, has no macroscopic effect on tiling entropy in the thermodynamic 
limit. The tilings become homogeneous, displaying the free-boundary entropy density throughout. 

The idea is to consider tilings (we focus on "diagonal" tilings with k\ = k 2 = fc 3 = p) that fill a regular octahedron 
O instead of the rhombic dodecahedron RD. Eight vertices of the RD must be truncated to produce the O. We call 
tilings with octahedron boundaries OS-tilings. Such an octahedron is displayed in figure 3. It is inscribed in an RD 
and has puckered boundaries instead of flat ones. Despite this puckering, the boundaries are effectively flat in the 
thermodynamic limit. The same kind of idea is developed in reference 14 in the hexagonal case where a puckered two- 
dimensional hexagonal boundary is introduced so that the tilings are homogeneous and exhibit no frozen regions. It is 
demonstrated (and numerically checked) that the random tilings filling this puckered hexagon display a free-boundary 
entropy in the thermodynamic limit. 
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FIG. 3. The puckered octahedral boundary conditions (fci = hi = £3 = p). The octahedron O is inscribed in the rhombic 
dodecahedron RD of side p coming from the solid partition method. The volume ratio |i?D|/|0| = 3/2. 

We now explain why OB-tilings have a free boundary entropy even though their boundary is fixed. The tiling 
entropy is contained in fluctuations of the associated hypersurface (p. Small-scale fluctuations are integrated in an 
entropy functional s[cj>] which takes into account the total number of possible piecewise linear hypersurfaces ip that 
are close to the smooth one <f> (see reference [ 14 ] for discussion). Functions maximizing s[<j>] represent the dominant 
macroscopic states of the system. Note that s[<j>] can be written as a functional of the gradients of <fi, known as the 
phason gradient or phason strain in the quasicrystal community. And s[(p] is maximum and equal to the free-boundary 
entropy (J free when this gradient vanishes everywhere. 

Fixed boundaries on tilings translate into fixed boundary conditions for the functions <j>. Therefore s[<f>] must be 
maximized on a restricted set of functions, F. For rhombic dodecahedral boundaries on RDB tilings, the boundaries 
of functions <f> £ F are non-flat polyhedra, and the phason gradient cannot vanish everywhere. Therefore their 
entropy density <Jfi xe d is bounded below by a free- For octahedral boundaries on OB tilings, the functions <j> are also 
constrained by a fixed boundary condition. But in this case, the boundary is flat and strain-free. It does not impose 
any phason gradient on the functions <j>. The phason gradient can vanish everywhere and the maximum of s[<f>] equals 

& free • 
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FIG. 4. Boundary conditions for the partition problem associated with the strain-free tiling problem with boundary conditions 
of Fig. 3. Two views of the truncated cube (of side p) are provided. The marked vertex is the same in both views. 

Figure 4 illustrates changes in the partition array boundary conditions needed to achieve OB tilings instead of RDB 
tilings. The partition array is no longer a cube, but two opposite pyramidal corners have been truncated, leaving a 
slab of D^d symmetry. This slab contains N p = p 3 — (p — l)p(p + l)/3 parts. 

We indicate minimal and maximal values of the parts: the three faces in the topmost figure bear minimal values of 
their adjacent boxes, whereas the three hidden faces bear maximal values. These values are visible on the bottom in 
the view in the bottommost figure. Individual parts are bounded by integers depending on the position of the part in 
the array. Since interior parts are bounded by their neighbors, it is sufficient to impose these bounds on surface parts, 
as indicated in figure 4. The lower bounds range from to p— 1, with constant values on lines parallel to diagonal of 
cube faces. The upper bounds are read on the opposite faces, and range from 1 to p. 

Using the methods discussed in section II A, such partitions generate tilings filling the octahedron O. For an original 
partition cube of sides k\ = &2 = ^3 = P, the full RD contains N t — 4p 3 tiles (see equation (3)). The octahedron 
O contains Nt = 4p 3 — 4(p — l)p(p + l)/3 tiles (these are the Ap 3 tiles in the RD minus the tiles in the truncated 
corners). These partitions are used in section III to estimate numerically o/ ree . 

C. Arctic octahedron conjecture (Linde, Moore and Nordahl) 

The relation between free and fixed boundary entropies in not known in general, because the generalization of the 
proof of the hexagonal case requires the knowledge of the free boundary entropy. In two dimensions, this relation is 
non-trivial because the unique function max maximizing s[(f>] has a complex expression 6 ' 14 , leading to the arctic circle 
phenomenon: the statistically dominating tilings are frozen (and periodic) outside a circle inscribed in the boundary 
hexagon. 

In three dimensions, Linde, Moore and Nordahl have recently explored numerically the typical shape of a RDB- 
tiling and have conjectured 15 that the two-dimensional circle then becomes not a sphere but a regular octahedron, 
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inscribed in RD like O in figure 3. More precisely, the unfrozen region is not exactly O but tends towards O at the 
large size limit. 

This conjecture has a crucial consequence 16 : the statistical ensemble of RDB-XAYmgs is dominated by tilings periodic 
outside O and random inside O, equivalently OS-tilings like in the previous section completed by eight periodically 
tiled pyramids to fill RD. Note that James Propp 17 has also conjectured independently of the present work that 
i?_D£>-tilings should be homogeneous inside O. Taking into account the tiles in the frozen regions to calculate an 
entropy per tile, one finally gets 

_ Np _2 

& fixed — AT - 0~free — ^"/ree> \^) 

since the ratio of the numbers of tiles in RD and O is 3/2. The next section is dedicated to the calculation of this 
ratio via Monte Carlo simulations. 



III. MONTE CARLO 



Conventional Monte Carlo simulations using Metropolis sampling 2 are useful for generating typical tiling configu- 
rations. They can reveal temperature driven phase transitions 25 and yield quantitative evaluations of phason elastic 
constants from fluctuations. Evaluation of entropy by Monte Carlo simulation demands specialized techniques. Ad- 
vanced methods using specialized dynamics based on entropic sampling 26 and transfer matrices 27 ' 28 yield reasonably 
accurate free energies and absolute entropies. 

We employ a variant of the transition matrix method 29 ' 30 that couples a conventional Metropolis Monte Carlo 
simulation with a novel data collection and analysis scheme to construct a numerical approximation to the transition 
matrix, described below in section III A. The density of states is an eigenvector of the exact transition matrix, and the 
sum of the density of states yields the total number of states. This method yields highly accurate absolute entropies 
with impressive efficiency. 



A. Transition Matrix 



For any legal partition P = {riijh}, wc define its "energy" as its total height 

E(P)=J2n ijk , (5) 

ijk 

and incorporate the constraint of legality by defining E = co for any partition P that violates the ordering conditions 
inside the partition or at the boundaries. Under single vertex flip dynamics (see Fig. 2) a single partition element 
increases or decreases by one unit resulting in an energy change AE = ±1. The ground state of this model is the 
lowest height legal partition. For the boundary conditions employed here, the ground state is unique and we denote 
its energy as E min . There is also a unique maximum energy state of energy E max . 

In a Metropolis Monte Carlo simulation of the partition, a randomly chosen partition element is randomly increased 
or decreased by one. The change is accepted if it lowers the energy (i.e. we attempted to decrease it, and the result 
is a legal partition). If the change raises the energy, it is accepted with probability exp (— E/T). Note that illegal 
partitions are never accepted. Low temperatures bias the partitions towards low heights, while high temperatures 
remove any bias according to height. As T becomes large, the Metropolis simulation faithfully reproduces the random 
tiling ensemble in which all configurations have equal weight. If we choose, we may take negative values of temperature 
T. Large negative temperatures again reproduce the random tiling ensemble, while small negative temperatures favor 
partitions of maximal height. 

For a given partition P, a number of partitions n±(P) can be reached by single upwards or downwards steps. 
Additionally, a certain number uq(P) of steps are forbidden due to the partition constraints. The sum rule 
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n-(P) + no(P) + n+(P) = 2N P 



(6) 



holds for every partition P. This value is twice the number of parts because each part can be raised or lowered. For 
any given partition P, the calculation of n± ; o is easy (though not fast) and exact. 

We define the transition matrix u)± t o(E) as a matrix with E max — E min + 1 rows (one for each allowed energy E) 
and three columns (labeled "+" for upwards transitions, for forbidden transitions and "-" for downwards transitions. 
Alternatively, we can think of these three columns as the diagonal and two off-diagonals of a square matrix u>(E, E') 
of dimension E max — E min + 1. Formally, we define 

u±,o(E) = E n ± , (P)/2N p (7) 

^ ' P(E) 

where the sum is over all partitions of energy E and the normalization W(E) is the total number of partitions with 
energy E. In general the set P(E) and the value W(E) are not known, preventing us from actually calculating the 
transition matrix. However, we may calculate the matrix numerically with high precision by averaging n±fi{P)/2N p 
over those partitions P occurring during a Monte Carlo simulation. By virtue of eq. (6), the matrix elements obey 
the sum rule 



u-(E) + u (E) + cj+(E) = 1 (8) 

for any energy E, so each row of the transition matrix can be normalized independently without knowledge of W(E). 

The transition matrix can be interpreted in terms of the rates at which Monte Carlo moves would be accepted in a 
hypothetical simulation at infinite temperature (even if our actual simulation is performed at finite temperature) . In a 
single Monte Carlo step, the probability for upwards or downwards transitions at infinite temperature is n±(P)/2N p , 
and the probability a move will be rejected is n${P) /2N p . We can predict finite temperature transition probabilities by 
multiplying u)± t o(E) with appropriate Boltzmann factors, and in fact these agree well with actual observed acceptance 
rates. 

To extract the density of states and the entropy, consider the infinite temperature detailed balance condition. At 
infinite temperature the probability a randomly chosen partition has energy E is W(E)/Z, where 

Z = Y J W{E) (9) 

E 

is the total number of legal partitions (the "partition function"). Multiplying the probability W(E)/Z by the accep- 
tance rate of upwards transitions lo + (E) yields the total forwards transition rate. The backwards rate is obtained in 
similar fashion. Detailed balance requires that the total rate of forward transitions from energy E to energy E + 1 
must equal the total rate of backwards transitions, hence 

lu + (E)W(E)/Z = u)-{E+ l)W{E + l)/Z. (10) 

It is useful to rearrange the detailed balance eq. (10) to find 

W(E + l) = "- {E + 1) W(E) (11) 
U)+{E) 

which allows us to iteratively extract the full density of states function W(E) using uniqueness of the ground state, 
W(E min ) = 1. Finally, the total entropy 

S = lnZ, (12) 

and the entropy density a — S/N t . 

Note that these infinite temperature transition rates may be calculated from finite temperature Metropolis simula- 
tions. Indeed, certain low energies E may occur very rarely in high temperature simulations so we must perform low 



8 



temperature simulations to generate partitions with energy E, from which the infinite temperature transition rates 
may be calculated. By sweeping over a range of temperatures we obtain accurate values of u>±(E) for all energies. 
Then, we use the infinite temperature detailed balance condition (10) to extract the density of states even though our 
simulation is performed entirely at finite temperatures. 

The accuracy of our result is controlled by the accuracy with which we determine w± (E) . For each partition, the 
transition numbers n±_o(P) are calculated exactly and added into the row of the transition matrix with energy E(P). 
Since there may be many partitions with the same energy E, each having different transition numbers, the accuracy 
with which a row of the transition matrix is determined is limited by our ability to generate a representative sample of 
partitions. We store the matrix as the integer valued sum of the integers n±,o(P), and impose the normalization (8) 
after data collection is complete. Because the values of n± t o(P) can be quite large, and we visit each energy a very 
large number of times, ordinary 4 byte integers can not hold the data. We implemented special procedures to handle 
storage and algebraic manipulations of large integers. 

B. Numerical Data 

Fixed boundary partitions are initialized at zero height (and thus E = E m i n = 0). Their maximum energy 
E m ax = k\kik-iV- We accumulate data in the transition matrix through four sweeps over temperature: from T min to 
T m ax during which time the mean energy grows from E min to about E max /2; from — T max to —T min during which 
time the mean energy grows from E max /2 to E max ; from — T m j„ back to —T max ; from T max back down to T m i n . Each 
sweep visits 201 temperatures in a geometric sequence from the initial to the final temperature. Free boundary tilings 
are initialized with a flat partition at energy close to (E min + E max )/2. We again perform four sweeps: from T m ax 
down to T min ] from T min back up to T rnax ; from -T max to -T min ; from -T min back to -T max . 

For both fixed and free boundary tilings, symmetries of the model dictate that the density of states is symmetric 
about the midpoint energy E mi d = {E min + E max )/2. The density of states W(E) has a strong maximum at the 
midpoint. By including both positive and negative temperature sweeps we sample both low and high energy states. 
By reversing our sweeps we mitigate possible systematic sampling errors associated with the direction of the sweep. 

At each temperature during a sweep, we accumulate data on Nmc sample configurations. For a given configuration 
(partition P), the time required to calculate n± ; o(-P) is proportional to the partition size N p , and hence to the 
time required to attempt N p flips. Because the data accumulation is time consuming, we can perform many single 
vertex flips between sample without significantly slowing down the simulation. We take Nfl — N p /5 single vertex 
flips, which yields rough equality between time spent flipping and collecting data. Prior to data collection at any 
temperature, we anneal at fixed temperature for Nmc x Nfl/100 steps. The total number of attempted flips in an 
Nmc = 10 6 run on a 10 x 10 x 10 partition is thus in excess of 4 x 201 x 10 6 x 10 3 /5 = 1.6 x 10 11 attempted flips. A 
run of this length takes 5 days on a 1.7GHz Pentium 4 processor. 

This protocol was chosen to ensure approximately uniform coverage of energies from E min to E max . The energy 
distribution of partitions visited at each temperature overlaps strongly the energy distribution of partitions visited 
at the previous and following temperatures. The value of T m i n is chosen sufficiently low to guarantee some coverage 
of the extreme energy states. Because there are few states at the energy extremes, there is an entropic barrier to 
reaching these extremes. This causes a spike in the coverage close to the extreme energies that cannot be avoided using 
Metropolis sampling. The value of T max is chosen sufficiently high that the midpoint energy is a local maximum in the 
coverage. Fig. 5 plots the number of partitions sampled as a function of energy during the longest runs (length 10 6 ) 
for the p = 4 partition. Note that the number of hits exceeds 10 6 uniformly for each energy in the range < E < 256. 

The density of states, W(E) shown in figure 5 (center), is nearly a Gaussian. W{E) reaches a peak value of 
1.7 x 10 15 states at energy E = 128, and has a half width at half maximum of AE — 20. As the system size grows 
this width grows more slowly than the number of tiles, so the density of states asymptotically approaches a delta 
function. Although we sampled a total of 8 x 10 6 configurations at E = 128, this represents a fraction of only about 
4.7 x 10~ 9 of the total number that exist at this energy. 

The microcanonical entropy, defined as logW(E), is plotted in the lower panel of Fig. 5. This plot reveals the 
expected symmetry around E mi d- The degree to which symmetry is broken can be used as an indicator of errors 
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accumulated during our iterative calculation (11), because sampling errors at low and high energies need not cancel. 
By inspection it can be seen that this error is quite low, since the entropy returns essentially to zero at high energy. 




Energy 

FIG. 5. Simulation data for p = 4 fixed boundary tiling. Top panel: Number of configurations sampled (units of 10 6 ). Center 
panel: Density of states W(E) (units of 10 14 per energy). Bottom panel: Microcanonical entropy S(E) = InW(E). 

Table I shows the convergence of entropy data as run length grows. We show the convergence only for the largest 
system size, which represents our worst case. The residual R = log W(E max ) measures the failure of the calculated 
microcanonical entropy to return to zero at high energy. It reflects the cumulative error in W(E), and thus, assuming 
statistical independence of the errors at each energy, it provides an upper bound on the error in W(E) for any energy. 
Converting the error in W(E) into an error in the entropy, we estimate an uncertainty 

|Aa| w R/N t . (13) 

Our numerical simulation results in Table I, and also simulations of smaller systems for which the entropy is known 
exactly, are in good numerical agreement with this estimate. 

Table II summarizes our data for all system sizes studied. We report the entropy resulting from the run of length 
Nmc = 10 6 , and use the difference between that run and the length Nmc — 10 5 run for the quoted uncertainty. 

To extrapolate our data to infinity, we fit to functional forms. For the fixed boundary entropy, we expect finite-size 
corrections of order \og(p) jp and 1/p because of boundary effects. Indeed, these two terms are the first correction terms 
in the exactly solvable one- and two-dimensional random tilings. Note that they have the same order of magnitude 
for the values of p considered here. This logarithmic form fits the data much better than a simple correction of order 
1/p. For fixed boundary tilings the data fit well to 

a fixed ( P )^ 0.145 -0.0049^ + ^, (14) 
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from which we conclude that a fixed — 0.145(3). The value is obtained from a fit excluding the data point with p = 1, 
while the uncertainty estimate comes from excluding instead the data point with p = 10. This value is very close to 
a conjectured limit 9 of 0.139, but given our small uncertainty, we believe the conjectured value is not exact. 
For the free boundary data we find 

a free (p)^ 0.214 - 0.052^-™^, (15) 

p p 

from which we conclude that cr/ ree = 0.214(2). As in the fixed boundary case, the logarithmic finite size correction 
fits the data better than a simple correction of order 1/p. Should the octahedral fixed boundary have a trivial 
short-range effect, it would create only a 1/p leading term (from surface vs bulk contributions), as indeed occurs in 
the corresponding two-dimensional hexagonal case 14 (see II B). Therefore, it appears that the strain- free octahedral 
boundary has a nontrivial long-range effect on the local entropy density, presumably by limiting the fluctations of the 
height function within a "penetration depth" that depends on the side length p. A penetration depth growing like 
log(p) is consistent with our leading finite-size correction. 
For the ratio, we fit to 

g *"<P> ~ 1.48 -0.37^-^, (16) 

VfixedKP) P P 

and we conclude the ratio a free/ & fixed — 1-48(3). This value equals 3/2 within our uncertainty. 

It would be desirable to achieve higher precision in this ratio in the future. At present, we are limited to system 
sizes of p = 10 or less because for larger systems the density of states W{E) exceeds 10 308 , the limiting floating point 
number on our Pentium 4 processor. Either specialized floating point arithmetic or a 64-bit processor will be required 
to treat systems of size p = 11 and above. 



IV. DISCUSSION 



Returning to the arctic octahedron conjecture, we recall (see equation (4)) we expected a ratio a free/ 'a fixed = 3/2. 
In fact, with few extra hypotheses, we now demonstrate that our numerical results provide a strong support to the 
arctic octahedron conjecture. Suppose that, as in the two-dimensional case 31 , the function 4> max maximizing the 
entropy functional s[<f>] subject to a given boundary condition is unique 14 . If cf>o is the piecewise linear function which 
vanishes identically in the central octahedral region O and has maximum strain (and entropy) in the eight pyramidal 
regions comprising RD — O, then s[<f> ] =2/3 a free-, and thus s[<f>o] = a fixed- Therefore <po = </> max , by uniqueness of 

0max ■ 

As a conclusion, if <J free/ & fixed = 3/2 and </> max is unique, then the arctic octahedron conjecture is true. Note that 
this conjecture can be generalized to higher dimensions, making in principle possible a calculation of the above ratio 
for any dimension, at least as far as codimension one problems are concerned. 

In terms of de Bruijn membranes, the above result means that all de Bruijn membranes are straight (or flat), at 
least at large scales. Indeed if de Bruijn membranes are flat in a i?D-tiling, the 4 de Bruijn families intersect in the 
central octahedron O, but only intersect 3 by 3 in the 8 pyramidal corners, leading to a frozen tiling in these 8 regions 
and a strain- free tiling in O. 

By comparison, onc-dimensional de Bruijn lines in two dimensional hexagonal tilings are not straight since that 
would lead to a hexagonal central arctic region with uniform entropy density. If crf^ denotes the free boundary 
entropy of hexagonal tilings, one would get a fixed boundary diagonal entropy per tile of 3 cr 3 ^/^ = 0.242 instead of 
0.261 (using similar arguments as above). Fluctuations of the de Bruijn lines raise the entropy from 0.242 to 0.261. 
Entropic repulsion of the lines causes them to bend. 

In contrast, the repulsion between de Bruijn membranes is sufficiently weak that they are not forced away from their 
flat configuration. This is possibly related to the fact that the fluctuation of a free 2-dimensional directed membrane 
in 3-dimensional space is of order ^/log(L) where L is its linear size 3 , whereas the fluctuation of a free 1-dimensional 
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directed polymer in 2-dimensional space is larger and of order \f~L. Therefore it is natural to suppose that the flatness 
of the arctic region will persist in dimensions higher than 3 where the fluctuations are even smaller, since they are 
bounded 3 . 

Our result emphasizes important dimensional dependence of the transition between the frozen and unfrozen regions. 
Indeed, in 2 dimensions, the transition is continuous, since the entropy density is by the arctic circle and then 
continuously varies to reach its maximum value near the center of the hexagon, with a non-zero gradient everywhere 
except near the center. By contrast, the situation seems to be radically new and different in 3 dimensions, since 
the entropy density appears to be constant in the arctic octahedron O, with a vanishing gradient everywhere and a 
discontinuous transition at the boundary of O. 

This result (as well as its possible generalization to higher dimensions) is a strong support in favor of the partition 
method, of which it was formerly believed that it could not easily provide relevant results about free boundary 
entropies. Indeed, provided the arctic region is polyhedral and its boundary is strain-free, the ratio of both entropies 
is nothing but a ratio of volumes of suitable polytopes. The latter two conditions should be fulfilled as soon as the 
entropic repulsion between de Bruijn membranes is sufficiently weak, that is as soon as the spatial dimension is 3 or 
greater. 

To finish with, we mention that this transition matrix Monte Carlo technique can easily be adapted to the numerical 
calculation of the entropy of two-dimensional rhombus tilings. Indeed, the structure of the configuration space of such 
tiling problems in terms of flips has been characterized in reference 12 . Note however that no such simple result as the 
"arctic octahedron phenomenon" is expected in these two-dimensional classes of tilings, but the calculation of their 
configurational entropy is a challenge in itself. This work is in progress. 
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TABLE I. Convergence of p = 10 entropy data for increasing run length. The residual R = log W(E ma x) measures the 
cumulative error in W{E). 



Nmc 


10 3 


10 4 


10 5 


10 6 


G fixed 


0.147827 


0.147385 


0.147400 


0.147349 


-Rfxixed 


3.474 


0.180 


0.429 


0.002 


@ free 


0.197516 


0.197300 


0.197234 


0.197273 


Rfree 


1.214 


0.143 


-0.202 


-0.036 



TABLE II. Size-dependent entropies. Values in parentheses are uncertainties in final digit. Values without uncertainties are 
exact . 



p 


& free 


G fixed 


& f reef & fixed 


1 


0.1732868 


0.1732868 


1.000 


2 


0.1732868 


0.1601239 


1.080 


3 


0.17947(2) 


0.1545769 


1.161 


4 


0.18455(6) 


0.1517949 


1.216 


5 


0.18829(3) 


0.15017(2) 


1.254 


6 


0.19108(4) 


0.14918(6) 


1.281 


7 


0.19320(4) 


0.14848(1) 


1.301 


8 


0.19486(2) 


0.14780(1) 


1.318 


9 


0.19618(1) 


0.14762(2) 


1.329 


10 


0.19727(4) 


0.14735(5) 


1.339 


oo 


0.214(2) 


0.145(3) 


1.48(3) 
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